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Gross-Oliveira-Kohn density functional theory (GOK-DFT) for ensembles is in principle very 
attractive, but has been hard to use in practice. A novel, practical model based on GOK-DFT for 
the calculation of electronic excitation energies is discussed. The new model relies on two mod¬ 
ifications of GOK-DFT: use of range separation and use of the slope of the linearly-interpolated 
ensemble energy, rather than orbital energies. The range-separated approach is appealing as it en¬ 
ables the rigorous formulation of a multi-determinant state-averaged DFT method. In the exact 
theory, the short-range density functional, that complements the long-range wavefunction-based en¬ 
semble energy contribution, should vary with the ensemble weights even when the density is held 
fixed. This weight dependence ensures that the range-separated ensemble energy varies linearly 
with the ensemble weights. When the (weight-independent) ground-state short-range exchange- 
correlation functional is used in this context, curvature appears thus leading to an approximate 
weight-dependent excitation energy. In order to obtain unambiguous approximate excitation ener¬ 
gies, we propose to interpolate linearly the ensemble energy between equiensembles. It is shown that 
such a linear interpolation method (LIM) can be rationalized and that it effectively introduces weight 
dependence effects. As proof of principle, LIM has been applied to He, Be, H2 in both equilibrium 
and stretched geometries as well as the stretched HeH"*" molecule. Very promising results have been 
obtained for both single (including charge transfer) and double excitations with spin-independent 
short-range local and semi-local functionals. Even at the Kohn-Sham ensemble DFT level, that is 
recovered when the range-separation parameter is set to zero, LIM performs better than standard 
time-dependent DFT. 


I. INTRODUCTION 

The standard approach for modeling excited states in 
the framework of density-functional theory (DFT) is the 
time-dependent (TD) linear response regime [1]. Despite 
its success, due to its low computational cost and rela¬ 
tively good accuracy, standard TD-DFT still suffers from 
various deficiencies, one of them being the absence of 
multiple excitations in the spectrum. This is directly 
connected with the so-called adiabatic approximation 
that consists in using a frequency-independent exchange- 
correlation kernel in the linear response equations. In or¬ 
der to overcome such limitations, the combination of TD- 
DFT with density-matrix- [2] or wavefunction-based [3-5] 
methods by means of range separation has been investi¬ 
gated recently. 

In this work, we propose to explore a time-independent 
range-separated DFT approach for excited states that is 
based on ensembles [6, 7]. One of the motivation is the 
need for cheaper (in terms of computational cost) yet 
still reliable (in terms of accuracy) alternatives to stan¬ 
dard second-order complete active space (CASPT2) [8] 


or N-electron valence state (NEVPT2) [9, 10] perturba¬ 
tion theories for modeling, for example, photochemical 
processes [11, 12]. Ensemble range-separated DFT was 
initially formulated by Pastorczak et al. [13] The authors 
considered the particular case of Boltzmann ensemble 
weights. The latter were controlled by an effective tem¬ 
perature that can be used as a tunable parameter, in ad¬ 
dition to the range-separation one. As shown in Ref. [14], 
an exact adiabatic connection formula can be derived for 
the complementary short-range exchange-correlation en¬ 
ergy of an ensemble. Exactly like in Kohn-Sham (KS) 
ensemble DFT [7, 15, 16], that is also referred to as 
Gross-Oliveira-Kohn DFT (GOK-DFT), the variation of 
the short-range exchange-correlation density functional 
with the ensemble weights plays a crucial role in the cal¬ 
culation of excitation energies [14]. So far, short-range 
density-functional approximations have been developed 
only for the ground state, not for ensembles. Conse¬ 
quently, an approximate (weight-independent) ground- 
state functional was used in Ref. [13]. 

The weight dependence of the range-separated ensem¬ 
ble energy and the ambiguity in the definition of an ap- 
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proximate excitation energy, that may become weight- 
dependent when approximate functionals are used, will 
be analyzed analytically and numerically in this work. 
By analogy with the fundamental gap problem [17], a 
simple and general linear interpolation method is pro¬ 
posed and interpreted for the purpose of defining un¬ 
ambiguously approximate weight-independent excitation 
energies. The method becomes exact if exact functionals 
and wavefunctions are used. The paper is organized as 
follows: After a brief introduction to ground-state range- 
separated DFT in Sec. II A, GOK-DFT is presented in 
Sec. IIB and its exact range-separated extension is for¬ 
mulated in Sec. IIC. The weight-independent density- 
functional approximation is then discussed in detail for a 
two-state ensemble. The linear interpolation method is 
introduced in Sec. IID and rationalized in Sec. HE. The 
particular case of an approximate range-separated en¬ 
semble energy that is quadratic in the ensemble weight 
is then treated in Sec. IIF. Comparison is made with 
Ref. [13] and time-dependent adiabatic linear response 
theory in Sec. IIG. A generalization to higher excita¬ 
tions is then given in Sec. IIH. After the computational 
details in Sec. Ill, results obtained for He, Be, H 2 and 
HeH+ are presented and discussed in Sec. IV. We con¬ 
clude this work with a summary in Sec. V. 


II. THEORY 

A. Range-separated density-functional theory for 
the ground state 

According to the Hohenberg-Kohn (HK) theorem [18], 
the exact ground-state energy of an electronic system can 
be obtained variationally as follows, 

Eq = min |F[n]+y dru„e(r)n(r)|, (1) 

where Une(r) is the nuclear potential and the minimiza¬ 
tion is performed over electron densities n(r) that inte¬ 
grate to a fixed number N of electrons. The universal 
Levy-Lieb (LL) functional [19] equals 

F[n] = min {'S\f + Wee\'S), (2) 

where f and IFee = kinetic en¬ 

ergy and regular two-electron repulsion operators, respec¬ 
tively. Following Savin [20], we consider the decomposi¬ 
tion of the latter into long- and short-range contributions, 

l/ri2 = w^;’^(ri2) -b<e’^(ri2), 

wiT{ri2) = erf(Airi2)/ri2, (3) 

where erf is the error function and /r is a parameter in 
[0,-|-oo[ that controls the range separation, thus leading 
to the partitioning 

F[n] = F^’^’>^[n]+E-;^^[n], (4) 


with 

= minl^-lf-bIFi^’^14'), (5) 

and IFee^ = The complementary /r- 

dependent short-range density-functional energy E^^ [n] 
can be decomposed into Hartree (H) and exchange- 
correlation (xc) terms, in analogy with conventional KS- 
DFT, 

^ y ydrdr'n(r)n(r'Xe’^(lr-r'l). (6) 

Inserting Eq. (4) into Eq. (I) leads to the exact expression 

Eo = mm {(TIT + + Kej'I') + i?H;^[n^]} 

= (4/^ If + + t4e 1 4 /^) + K.], (7) 

where Fie = / dr Vne(r) h(r) and h(r) is the density oper¬ 
ator. The electron density obtained from the trial wave- 
function 4t is denoted nip. The exact minimizing wave- 
function 4 'q in Eq. (7) has the same density as the 
physical fully-interacting ground-state wavefunction 4^0 
and it fulfils the following self-consistent equation: 

= ( 8 ) 

where 

/ X r 1 

dr h(r). (9) 

dn[r) 

It is readily seen from Eqs. (3) and (8) that the KS and 
Schrodinger equations are recovered in the limit of /i = 0 
and fj, —>■ -boo, respectively. An exact combination of 
wavefunction theory with KS-DFT is obtained in the 
range of 0 < /x < -boo. 

In order to perform practical range-separated DFT calcu¬ 
lations, local and semi-local short-range density function¬ 
als have been developed in recent years [21-24]. In ad¬ 
dition, various wavefunction-theory-based methods have 
been adapted to this context in order to describe the long- 
range interaction: Hartree-Fock (HF) [25, 26], second- 
order Mpller-Plesset (MP2) [25, 27, 28], the random- 
phase approximation (RPA) [29, 30], configuration in¬ 
teraction (GI) [31, 32], coupled-cluster (GG) [23], the 
multi-configurational self-consistent field (MGSGF) [26], 
NEVPT2 [33], one-electron reduced density-matrix- 
functional theory [34] (RDMFT) and the density ma¬ 
trix renormalization group method [35] (DMRG). In this 
work, GI will be used. The orbitals, referred to as HF 
short-range DFT (HF-srDFT) orbitals in the following, 
are generated by restricting the minimization on the first 
line of Eq. (7) to single determinantal wavefunctions. 
Note that, when /x = 0, the HF-srDFT orbitals reduce 
to the conventional KS ones. 

Finally, in connection with the description of excited 
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states, let us mention that the exact auxiliary excited 
states {^'f}*>o that fulfil the eigenvalue equation, 

( 10 ) 

can be used as starting points for reaching the physical 
excitation energies by means of extrapolation tech¬ 
niques [36-38], perturbation theory [39], time-dependent 
linear response theory [4, 5] or ensemble range-separated 
DFT [13, 14], as discussed further in the following. 


B. Ensemble density-functional theory for excited 
states 

According to the GOK variational principle [6], that 
generalizes the seminal work of Theophilou [40] on 
equiensembles, the following inequality 

< Tr t^H , (11) 

where H = T + Wee + Ine and Tr denotes the trace, 
is fulfilled for any ensemble characterized by a set of 
weights w = (wo, wi, ..., wm-i) with wq > wi > ... > 
wm -1 > 0 and a set of M orthonormal trial wavefunc- 
tions {'I'fc}o<fe<M-i from which a trial density matrix 
can be constructed: 

M-l 

(12) 

k^O 

The lower bound in Eq. (11) is the exact ensemble energy 

M-l M-l 

= E = E (13) 

k=0 fc=0 

where is the exact fcth eigenfunction of H and Eq < 
El < ... < Em-i- In the following, the ensemble will al¬ 
ways contain complete sets of degenerate states (referred 
to as ’’multiplets” in Ref. [7]). An important consequence 
of the GOK principle is that the HK theorem can be ex¬ 
tended to ensembles of ground and excited states [7], thus 
leading to the exact variational expression for the ensem¬ 
ble energy, 

E'^ = min |F'^[n] + J dr z;ne(r)»^(r)|’ , (14) 

where the universal LL ensemble functional is defined as 
follows, 

F"'[n] = min (Xr \t'^{f + II4e)j | ■ (15) 

fw_>„ L L J J 

The minimization in Eq. (15) is restricted to ensemble 
density matrices with the ensemble density n: 

Tr f'^h(r) = np„(r) = n(r). (16) 


Note that, in the following, we will use the convention 
Wfc = 1 so that the ensemble density integrates 
to the number of electrons N. The minimizing density 
in Eq. (14) is the exact ensemble density of the physical 
system n"'(r) = 


In standard ensemble DFT [7], that is referred to as 
GOK-DFT in the following, the KS partitioning of the 
LL functional is used, 

E'-[n]=Tr[n]+E^^M, ( 17 ) 

where the non-interacting ensemble kinetic energy is de¬ 
fined as 

T^[n] = min | Tr |, (18) 

fw_>„ 1 L J J 

and E^^^[n] is the w-dependent Hxc functional for the 
ensemble, thus leading to the exact ensemble energy ex¬ 
pression, according to Eq. (14), 

E^ = inm {Tr + Ke)] + . (19) 

The minimizing GOK density matrix, 

M-l 

fr= ( 20 ) 

k=0 

reproduces the exact ensemble density of the physical 
system, 

np„(r) = n'^(r), (21) 

and it fulfils the stationarity condition (i£’^[r^] = 0 

where 

^W[pW] ^ ^ ^ 

M-l 

+ ^ u;fc£:r(l-(^fcl^fe)). (22) 

k=0 

The coefficients are Lagrange multipliers associated 
with the normalization of the trial wavefunctions 4*^ from 
which the density matrix is built. Gonsidering variations 
'kfe — t 'kfc -I- (5'1'fe for each individual states separately 
leads to the self-consistent GOK equations [7]: 




= £:ri$n, 0<fc<M-l. (23) 


C. Range-separated ensemble density-functional 
theory 

In analogy with ground-state range-separated DFT, 
the LL ensemble functional in Eq. (15) can be range- 
separated as follows [13, 14], 

E'" [n] = [n] + [n], (24) 
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where 


F 


= min {Tr {f + } . (25) 


In the following, the short-range ensemble functional 
will be partitioned into w-independent Hartree and w- 
dependent exchange-correlation terms, 

E^H.r[n]=E-’^[n]+E-’>^’-[n]. (26) 

Note that the decomposition is arbitrary and can be 
exact or not, depending on the short-range exchange- 
correlation functional used. In practical calculations, lo¬ 
cal and semi-local exchange-correlation functionals may 
not remove the so-called ’’ghost interactions” [41, 42] 
that are included into the short-range Hartree term. 
Such interactions are fictitious and unwanted. Their de¬ 
tailed analysis, in the context of range-separated ensem¬ 
ble DFT, is currently in progress and will be presented 
in a separate work. 

Combining Eq. (14) with Eq. (24) leads to the exact 
range-separated ensemble energy expression 

f-(f+ IT^^’^ + Ke) 

(27) 

The minimizing long-range-interacting ensemble density 
matrix f^’'" = reproduces the 

physical ensemble density. 


E^ = rninj Tr 

f'w L 


or as the slope of the linear interpolation between w = 0 
and w = 1/2, 


uj = 2(E^=^''^-Eo). (33) 


Let us stress that Eqs. (32) and (33) are equivalent in the 
exact theory. By using the decomposition (see Eqs. (27) 
and (28)) 

E^ = (l- w)(^')’“|f + 

+a;(4'f“|T + (34) 

that can be rewritten in terms of the auxiliary long-range 
interacting energies as follows, according to Eq. (29), 

E™ = (1 - 

- J dr ^ n-(r) + (35) 

where the physical ensemble density equals the auxiliary 
one (see Eq. (28)), 


n“'(r) = (1 - w)n^M.™(r) -|- w (r), (36) 


and by applying the Hellmann-Feynman theorem, 

|n,i,M.™(r), (37) 


dSr d 


^Hxc 

5n(r) 


we finally recover from Eq. (32) the following expression 
for the first excitation energy [14], 


nrM,w(r) = n'^(r), (28) 

and, by analogy with Eq. (22), we conclude that it should 
fulfill the self-consistent equation 

+ Ke + I [4/^’-) 

0<fc<M-l. (29) 

Note that the Schrodinger and GOK-DET equations are 
recovered for fj, —>■ -too and /r = 0, respectively. 


In the rest of this work we will mainly focus on ensem¬ 
bles consisting of two non-degenerate states. In this case, 
the ensemble weights are simply equal to 

wi = w, Wo = 1 — w, (30) 

where 0 < w < 1/2, and the exact ensemble energy is a 
linear function of w, 


= (l-w)Eo + wEi. (31) 


Consequently, the first excitation energy uj = Ei — Eq 
can be written either as a first-order derivative, 


dE“ 

U) = 


(32) 




dEl 


ST,fJ.,W [■ 

Hxc I- 
dw 


= A£f^’'' 


^xc • 


(38) 


It is readily seen from Eq. (38) that the auxiliary ex¬ 
citation energy differs in princi¬ 

ple from the physical one. They become equal when 
fjL —)■ - 1 - 00 . Eor finite /i values, the difference is simply 
expressed in terms of a derivative with respect to the en¬ 
semble weight = 9E®(.’^’’"[n]/9wl„^„„. Note that 
the Hartree term does not contribute to the second term 
on the right-hand side of Eq. (38) since it is, for a given 
density n, w-independent (see Eq. (26)). Interestingly, 
when w —)■ 0, an exact expression for the physical excita¬ 
tion energy is obtained in terms of the auxiliary one that 
is associated with the ground-state density (see Eq. (10)), 


LjJ 


= + 




dw 


(39) 


w—0 


Note also that, when ^ = 0 and the first excitation is a 
one-particle-one-hole excitation (single excitation), the 
GOK expression [7] is recovered from Eq. (38), 


w = Ae™ -f A 


W 

XC’ 


(40) 


where Ae“ is the HOMO-LUMO gap for the 

non-interacting ensemble and A]/. = dE^^[n]/dw\^^^„ ■ 
In the w —>■ 0 limit, the exact excitation energy can be 


dw 
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expressed in terms of the KS HOMO eo and LUMO Si 
energies as follows, 

UJ = - £0, (41) 

where = Si + As shown analytically by 

Levy [43] and numerically by Yang et al. [15], corre¬ 
sponds to the jump in the exchange-correlation potential 
when moving from w = 0 (ground state) to w > 0 (en¬ 
semble of ground and excited states). This is known as 
the derivative discontinuity (DD) and should not be con¬ 
fused with the ground-state DD that is related to ion¬ 
ization energies and electron affinities, although there 
are distinct similarities at a formal level [44-46]. Con¬ 
sequently, the quantity A[(^“ introduced in Eq. (38) will 
be referred to in the following as short-range DD. 


readily seen from Eq. (45). In other words, the short- 
range DD is not modeled at the WIDFA level of approxi¬ 
mation. Finally, the exact (/i-independent) ground-state 
energy will still be recovered when w —)■ 0 if no ap¬ 
proximation is introduced in the short-range exchange- 
correlation functional, 

= Eg. (46) 

Obviously, the exact ensemble energy will in general not 
be recovered for w > 0. By rewriting the WIDFA ensem¬ 
ble energy as 

E>^’^ = (1 - + raff 

dr (47) 

and applying the Hellmann-Feynman theorem. 


D. Weight-independent density-functional 
approximation and the linear interpolation method 


d£t 



\ 5n(r) J 


n^M.™(r), 


(48) 


Even though an exact adiabatic-connection-based ex¬ 
pression exists for the short-range ensemble exchange- 
correlation functional (see Eq. (133) in Ref. [14]), it 
has not been used yet for developing weight-dependent 
density-functional approximations. Let us stress that 
this is still a challenge also in the context of GOK- 
DFT [15]. A crude approximation simply consists in us¬ 
ing the ground-state functional [13], 

E-’'^’-[n]^£;->^[n], (42) 

thus leading to the approximate ensemble energy expres¬ 
sion 


E/^’^ = (1 - 

(43) 

that may depend on both /i and w, and where the ap¬ 
proximate auxiliary ensemble density equals 

h''’“’(r) = (1 - ■u;)n|,j..™ (r) +wnq,i.,^{r), (44) 

with 

i7''[h^’“']l^f’“) i = 0,l. (45) 

In the following we refer to this approximation as weight- 
independent density-functional approximation (WIDFA). 
Note that, at the WIDFA level, the ground-state density- 
functional Hamiltonian H^^[n] (see Eq. (9)) is used. The 
auxiliary wavefunctions associated with the bi¬ 

ensemble (0 < ic < 1/2) will therefore deviate from 
their ’’ground-state” limits 4'/ {w = 0) because of the 
ensemble density that is inserted into the short- 

range Hxc potential. Note that Eq. (45) should be 
solved self-consistently. Let us also stress that the 
ground-state short-range Hxc density-functional poten¬ 
tial i5A/[^(([n°]/(5n(r) is recovered in the limit w —)■ 0, as 


we see that, within WIDFA, the first-order derivative of 
the ensemble energy reduces to the auxiliary excitation 
energy that is in principle rc-dependent, 

AEn,w _ _ 

= Af. (49) 

Therefore, in practical calculations, the WIDFA ensem¬ 
ble energy may not be strictly linear in w, as illustrated 
for He in Fig. 1. In the same spirit as Ref. [17], we pro¬ 
pose to restore the linearity by means of a simple linear 
interpolation between the ground state (w = 0) and the 
equiensemble [w = 1/2), 

F'"’™ =EgE 2w(F''4/2 _ (50) 

This approach, that will be rationalized in Sec. H E, is 
referred to as linear interpolation method (LIM) in the 
following. The approximate excitation energy is then un¬ 
ambiguously defined as 


(jj 


u _ 

LIM ~ 


dF^’ 

dru 


2(F^4/2 _ ^g)_ 


(51) 


Note that, according to Eq. (33), LIM becomes exact 
when the exact weight-dependent short-range exchange- 
correlation functional is used. By analogy with the grand 
canonical ensemble [17], we can connect the linear inter¬ 
polated and curved WIDFA ensemble energies as follows, 

nW 

^ + df A^g^ (52) 

^0 

so that, according to Eqs. (49) and (51), 

= + (53) 

As readily seen from Eqs. (38) and (53), A^y plays the 
role of an effective DD that corrects for the curvature of 
the WIDFA ensemble energy, thus ensuring strict linear¬ 
ity in w. A graphical representation of LIM is given in 
Fig. 2. 
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E. Rationale for LIM and the effective DD 

The effective DD has been introduced in Eq. (52) 
for the purpose of recovering an approximate range- 
separated ensemble energy that is strictly linear in w. 
This choice can be rationalized when using a range- 
dependent generalized adiabatic connection formalism 
for ensembles (GAGE) [14], where the exact short-range 
ensemble potential is adjusted so that the auxiliary en¬ 
semble density equals the (weight-independent) density 
n(r) for any weight ^ and range-separation parameter ly 
values: 


the linearly-interpolated WIDFA ensemble energy 
is, in principle, as relevant as other choices. Still, the ana¬ 
lytical derivations and numerical results presented in the 
following suggest that LIM has many advantages from a 
practical point of view. By doing so, we hnally recover 
the expression in Eq. (53): 

(61) 

F. Effective DD and excitation energy for a 
quadratic range-separated ensemble energy 


+ y drv‘^’«(r)h(r)^ 

= £:r«K’«), * = 0,1, (54) 

where 

(1 - (r) + W;--? W = (55) 

It was shown [14] that the exact short-range ensem¬ 
ble exchange-correlation density-functional energy can be 
formally connected with its ground-state limit (w = 0) 
as follows, 

pW 

[n] = [n] + / d^ [n], (56) 

where the exact density-functional DD equals 

[n] = ■ (57) 

When rewritting the WIDFA ensemble energy in Eq. (43) 
as 

-f J dr v„e(r)h''’“’(r), (58) 

it becomes clear, from Eqs. (52) and (56), that LIM im¬ 
plicitly defines an approximate weight-dependent short- 
range exchange-correlation functional: 



In order to connect the exact DD with the effective one, 
let us consider Eq. (57) in the particular case n = 
and ^ = w, thus leading to 

ABr.M.»[^M.»] = A£+“’“’ - (60) 

where A£+°°’“’ is the excitation energy of the fully- 
interacting system with ensemble density h^’™. If the 
latter is a good approximation to the true physical ensem¬ 
ble density n™, which is the basic assumption in WIDFA, 
then becomes w-independent and equals the 

true physical excitation energy. As discussed previously, 
the latter has various approximate expressions that all 
rely on various exact expressions. Ghoosing the slope of 


For analysis purposes we will approximate the WIDFA 
ensemble energy by its Taylor expansion through second 
order in w (around w = 0) over the interval [0,1/2], 


= Eo+wE^^^^'> -f 


where, according to Eqs. (10), (45), (48) and (49), 






dw 


_ CM cM 

— C-y Cq , 


(62) 


(63) 


W—Q 


and 


^^(2) _ 




drdr' 




vj—0 

2 ZT-sr,^!- 01 


(5n(r')(5n(r) 


{n^^{r) -n°{r)) 


X n^^(r') - n°(r') -|- 


dnq,^.u,{E) 


dw 


w—0 


(64) 


As shown in Sec. IV, this approximation is accurate when 
^ > I.Ooq^. For smaller values, and especially in the 
GOK-DFT limit (/* = 0), the WIDFA ensemble energy 
is usually not quadratic in w. Nevertheless, making such 
an approximation gives further insight into the LIM ap¬ 
proach, as shown in the following. From the equiensemble 
energy expression 

^m.1/2 ^ (65) 

2 8 

and Eq. (51), we obtain the LIM excitation energy within 
the quadratic approximation, that we shall refer to as 
LIM2, 


<IM2 = - Ao) 


= ^m(i) _ 


( 66 ) 


thus leading to 


OJ- 


LIM2 

1 

'i 


= 


drdr 




Sn{r')Sn{r) 




X ( n^('(r') - n°(r') -|- 


5n,jM.™(r') 


dw 


W—0 


(67) 
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As shown in Appendix A, an explicit expression for 
the linear response of the ground-state density to 

variations in the ensemble weight w can be obtained from 
self-consistent perturbation theory. Thus we obtain the 
following expansion through second order in the short- 
range kernel: 

, _ CM _ pM 

^LIM2 “ <^1 “"0 

+ I / / drdr' ffHff| (n^.(r')-n°(r')) 

A J J 6n{r')6n(j) ^ ^> 

x(n,„ 5 *(r) -n°(r)) 


where 


dridr^drdr' 


(5n(r'j^)dn(ri) 


X2 psr.^tr 01 

xK.«)-nV.))E=^#d|(tl 


The latter expression is convenient for comparing LIM 
with time-dependent range-separated DFT, as discussed 
further in the following. Returning to the quadratic en¬ 
semble energy in Eq. (62), its Hrst-order derivative equals 

^ (69) 

aw 

thus leading to the following expression for the effective 
DD, according to Eq. (66), 


^eff — 


dE^’“ 

dw 


= (70) 

In conclusion, the effective DD is expected to vanish at 
w = 1/4 when the WIDFA ensemble energy is strictly 
quadratic, as illustrated in Fig. 2. 


AE(o) = - si; + E-CK.] - 

-f J h ri°{r) -n^^ir)), (74) 


and, according to Eq. (48), 

dA^ .fdtf: 

dw 




(5n(r) 

5n{r) 

(r) 



it is readily seen that the excitation energy will vary lin¬ 
early with w in the vicinity of w = 0. Therefore, in 
practical calculations, an optimal value for w must be 
determined [13]. This scheme can be compared with 
LIM2 by expanding the excitation energy in the density 
difference nq,i^{r) — n°(r), thus leading to 

AE{w) = - s; 

x(nvi,M(r) -n°(r)) 

9n|,^.£(r) 

+-. (7d) 

«=o 

or, equivalently, 

AE{w) = - s; 

+ l f f drdr' ^"fH^^°| ^.(r') - n°(r')) 

A J J 6n{r')6n(r) ^ 

.J^ N , 9 h^’“’«(r)| \ 


X (r) - (r)-I- 


G. Comparison with existing methods 

1. Excitation energies from individual densities 


where 


^A7,™.?(r) = (4u;-h^)n*M.«(r) -^n*M.£(r). (78) 


Pastorczak et al. [13] recently proposed to compute 
excitation energies as differences of total energies, 

AE{w) = Ei{w) - Eo{w), (71) 

where the energy associated with the state i (i = 0,1) is 
obtained from its (individual) density as follows: 

E,{w) = 

(72) 

From the Taylor expansion 

AE{w) = AE(0) + w +(!1(u;2), (73) 


This expression is recovered from the LIM2 excitation 
energy in Eq. (67) by applying the following substitution: 

n,j,g,j(r)h'^’“’^(r). (79) 

In other words, for a given ensemble weight w, the re¬ 
sponse of is used rather than the ground-state den¬ 

sity response in the calculation of the excitation energy 
AE{w). Note that integrating over space gives 

AwN. Therefore, may be considered as a density 

only when u; = 1/4. In this case, it is simply expressed 
as 

nA*d/4.?(j.) ^ (1 _^„-^_5(r), (80) 



and its response to changes in ^ equals 


and the gradient density vector becomes 






i=o 


= n^M(r) - n°(r) + 


dn 






(81) 


«=o 


Consequently, the LIM2 excitation energy can be recov¬ 
ered only if around ^ = 0, that means when 

the excitation energy reduces to the auxiliary one. Note 
finally that the averaged density in Eq. (80) can be in¬ 
terpreted as an ensemble density only if— 

It is unclear if its derivative at ^ = 0 has any physical 
meaning. 




' <*(!•) ' 

-<(r) 


( 88 ) 


The transition matrix elements associated with the 
density operator Ugj(r) have already been introduced in 
Eq. (A8). 


We propose to solve Eq. (82) by means of perturba¬ 
tion theory in order to make a comparison with LIM2. 
The perturbation will be the short-range kernel. Let us 
consider the auxiliary linear response equation, 


2. Time-dependent adiabatie linear response theory 


(4"]^ + - cc(a) X{a) = 0, (89) 


An approximation w to the hrst excitation energy can 
also be determined from range-separated DFT within the 
adiabatic time-dependent linear response regime [4, 5]. 
The associated linear response vector X fulhls 

-b X = 0, (82) 

where the long-range interacting Hessian and the metric 
equal 


Tjq — 




and 


S'[2]m = 


[Ri, ffjo 


[Ri, R 


■{[R^, 


R 


■j\o 


i\o 


-{\R^,R}]o 


,(83) 


(84) 


respectively. Short-hand notations [A, B]o = 
(T^|[i,H]|4/^), iLO = H^[n% and r] = |vl/f)(vl/^| 
with f > 0 have been used. The short-range kernel 
matrix in Eq. (82) is written as 


that reduces to Eq. (82) in the a = 1 limit, and the 
perturbation expansions 

X{a) = -b -b 0{a^), 

uj{a) = -b -b -b 0{a^). (90) 

Since we are here interested in the first excitation energy 
only, we have 

ri 


[0. 

Inserting Eq. (90) into Eq. (89) leads to the following 
excitation energy corrections through second order, 

^(1) 

w(2) =X(o)tx^^^^X(i), (92) 

where the intermediate normalization condition 
= 1 has been used, and 


, (91) 


jy-ST,n _ 


dr dr' 


,d^E\ 




Hxc 




Sn{r')Sn{r 
where the gradient density vector equals 


nW''(r')nW'^t(i.)^(85) 


j[i]/^(r) = 


[.R*,h(r)]o 


( 86 ) 


Since we use in this section a complete basis of orthonor¬ 
mal X-electron eigenfunctions {4'^}fc^o,i.... associated 
with the unperturbed long-range interacting Hamiltonian 
H>^[n^] and the energies }fe=o,i,..., orbital rotations 
do not need to be considered, in constrast to the ap¬ 
proximate multi-determinant formulations presented in 
Refs. [4, 5], such that matrices simply reduce to 


E, 


[An 




0 ■ 

0 ’ 


(87) 


-ba;(i)S'[2l^X(°). (93) 


According to Eqs. (85), (88) and (91), the first-order cor¬ 
rections to the excitation energy and the linear response 
vector become 




drdr 




n\ Q 

HxcL"- J /r.r\„n 


Sn{r')6n{r) 


<i(i’Xi(r), (94) 


and 


xw = - 


drdr' 




Sn{r')Sn{r) 




X ^ (n[il^(r') - <i(r')X(°^) j(95) 


respectively. Combining Eq. (85) with Eqs. (92) and (95) 
leads to the following expression for the second-order cor- 
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rection to the excitation energy: 


,( 2 ) = 


I I I I 



dn{r')6n{r) 


7iSi-.Mr„0l 


^ (ri)ng,(r ') 

s^-st 


+ E 


<^(I^lX^(I•') \ 
2Sl^-St-S^)- 


(96) 


with 


0 < w < , 

“ ~ Ml 




L=0 


( 102 ) 


and Eq < El <...< Ej are the 1 + 1 lowest energies 
with degeneracies {gL}o<L<i- In the exact theory, the 
ensemble energy is linear in w with slope 


The second summation in Eq. (96) is related to de¬ 
excitations. Within the Tamm-Dancoff approximation 
the latter will be dropped, thus leading to the following 
expansion through second order, according to Eqs. (91) 
and (94), 


dw 


= giEi 


Mi_i 


2 ^ gKJ^K , 
K=0 / 


thus leading to the following expression for the exact Jth 
excitation energy 


w = 


drdr 


dridr^drdr' 


Sn{r')6n{r) 

^Hxc[ 

5n(rj)(5n(ri) 


XCr'XiW 


A2psr,Mr„0l 

• ^ -^HxcF i M {r-)n^ fr' 1 

■5n(rO<In(r)”oiWnoi(rij^ 


X(riX»(j^') 


(97) 


A direct comparison can then be made with the LIM2 
excitation energy in Eq. (68). Thus we conclude that 
LIM2 can be recovered through first and second orders 
in the short-range kernel from adiabatic time-dependent 
range-separated DFT by applying, within the Tamm- 
Dancoff approximation, the following substitutions. 


''01 




(98) 


and 


E 

i>l 





^ 2 


X(ri)X(r') 


i>l 


_ 0i\ 


(99) 


respectively. 


LOj = Ej — Eq 


1 dEf 
gi dw 


1 

Mi-i 


i-i 

gK^K- 

K=1 


(104) 


The LIM excitation energy, that has been introduced 
in Eq. (51) for non-degenerate ground and first-excited 
states, can therefore be generalized by substituting the 
approximate first-order derivative (that may be both ji- 
and w-dependent) with its linear-interpolated value over 
the segment [0,1/M/], 


dEX 

dw 


Mi 


, (105) 


so that the Jth LIM excitation energy can be defined as 
Mi 


— 

Aim,/ 


gi 


(X’ 


IjMj — I 

— ^l-i 


1 


7-1 




(106) 


K=l 


where the equality = 7/’° has been used. In 

other words, LIM simply consists in interpolating linearly 
the ensemble energy between equiensembles that are de¬ 
scribed at the WIDFA level of approximation. 


H. Generalization to higher excitations 


III. COMPUTATIONAL DETAILS 


Following Gross et al. [7], we introduce the generalized 
w-dependent ensemble energy 


EY = X 

" M/_i 


gKEKj EwgjEj, ( 100 ) 


, K=0 


that is associated with the following ensemble weights, 

Wk= I (101) 

[w M/_i <k< Mj — 1, 


Eqs. (45) and (51) as well as their generalizations 
to any ensemble of ground- and excited states (see 
Eq. (106)) have been implemented in a development ver¬ 
sion of the DALTON program package [47, 48]. For sim¬ 
plicity, we considered spin-projected (singlet) ensembles 
only. In the latter case, the GOK variational principle is 
simply formulated in the space of singlet states [15]. In 
practice, both singlet and triplet states have been com¬ 
puted but, for the latter (that can be identified easily in a 
GI calculation), the ensemble weight has been set to zero. 
Both spin-independent short-range local density [20, 21] 
(srLDA) and Perdew-Burke-Ernzerhof-type [23] (srPBE) 
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approximations have been used. Basis sets are aug-cc- 
pVQZ [49, 50]. Orbitals relaxation and long-range corre¬ 
lation effects have been treated self-consistently at the 
full Cl level (FCI) in the basis of the (ground-state) 
HF-srDFT orbitals. For Be, the Is orbitals were kept 
inactive. Indeed, in the standard wavefunction limit 
{fj, —)■ -boo), deviations from time-dependent CC with 
singles and doubles (TD-CCSD) excitation energies are 
0.4 and 2.0 mEh for the 2s —)■ 3s and (2s)^ —)■ (2p)^ 
excitations, respectively. Comparisons are made with 
standard TD-DFT using LDA [51], PBE [52] and the 
Coulomb attenuated Becke three-parameter Lee-Yang- 
Parr [53](CAM-B3LYP) functionals. We investigated 
the following ensembles consisting of two singlet states: 
{1^5', 2^5'} for He and Be, {1^I]+, 2^E+} for the stretched 
HeH+ molecule and {1^11+, 2^S+} for H 2 at equilibrium 
and stretched geometries. For Be, the four-state ensem¬ 
ble {l^S”, 2^5', l^H} in Ag symmetry {l^D is doubly de¬ 
generate) has also been considered in order to compute 
the l^S' —> 1^1? excitation energy. 

IV. RESULTS AND DISCUSSION 
A. Effective derivative discontinuities 

1. GOK-DFT results (p = Q) for He 

Let us first focus on the GOK-LDA results (^ = 0 
limit) obtained for He. As shown in the top left-hand 
panel of Fig. 3, the variation of the auxiliary excitation 
energy with w is very similar to the one obtained at the 
quasi-LDA (qLDA) level by Yang et al. (see Fig. 11 
in Ref. [15]). An interesting feature, observed with 
both methods, is the minimum around w = 0.01. The 
derivation of the first-order derivative for the auxiliary 
excitation energy is presented in Appendix B. As readily 
seen from the expression in Eq. (BIO), at ru = 0, the 
derivative contains two terms. The first one, that is 
linear in the Hxc kernel, is expected to be positive 
due to the Hartree contribution. The second one is 
quadratic in the Hxc kernel and is negative (because of 
the denominator), exactly like conventional second-order 
contributions to the ground-state energy in many-body 
perturbation theory. The latter term might be large 
enough at rc = 0 so that the auxiliary excitation energy 
decreases with increasing w. The linearity in w (last 
term on the right-hand side of Eq. (BIO)) explains why 
that derivative becomes zero and is then positive for 
larger w values. As the excitation energy increases, the 
denominator mentioned previously also increases. The 
derivative will therefore increase, thus leading to the 
positive curvature observed for the auxiliary excitation 
energy. All these features are essentially driven by the 
response of the auxiliary excited state to changes in 
the ensemble weight (not shown). Returning to the top 
panels in Fig. 3, we see that the minimum at w = 0.01 
only appears when auxiliary energies are computed 


self-consistently. This is consistent with Eq. (BIO) 
where the second (negative) term on the right-hand side 
describes the response of the KS orbitals to changes in 
the Hxc potential through the w-dependent ensemble 
density. When the latter term is neglected, the auxiliary 
excitation energy has positive slope already at re = 0. 
For larger w values, self-consistency effects on the slope 
are reduced. Indeed, the response of the GOK orbitals is 
expected to be smaller as the auxiliary excitation energy 
increases. The large deviation of the non-self-consistent 
auxiliary excitation energy from the self-consistent one 
is due to the fact that, for the former, the ensemble 
density is constructed from the ground-state KS or¬ 
bitals. Finally, we note that the self-consistent auxiliary 
excitation energy equals the reference FCI one around 
w = 0.4. A very similar result has been obtained at the 
qLDA level by Yang et al. [15] We also find that both 
LDA and PBE yield very similar results. 

Let us now turn to the LIM excitation energy for 
/i = 0. By construction, it is w-independent, like in the 
exact theory. Note that the auxiliary excitation energy 
equals the LIM one for a w value that is slightly larger 
than 1/4, thus showing that the ensemble energy is not 
strictly quadratic in w. Moreover, as expected from the 
analysis in Appendix C, the effect of self-consistency is 
much stronger on the auxiliary excitation energy than 
on the LIM one. For the latter it is actually negligible. 
Turning to the effective DDs in the top panels of Fig. 3, 
these qualitatively vary with the ensemble weight similar 
to the accurate DD shown in Fig. 7 of Ref. [15]. Still, 
there are significant differences. For re = 0, the effective 
DD equals 0.0736 and 0.0814 at the LDA and PBE 
levels, respectively. The accurate value obtained by Yang 
et al. [15] is much smaller (0.0116 Eh). In addition, both 
LDA and PBE effective DDs equal zero close to w = 1/4 
that is much smaller than the accurate value of Ref. [15] 
(w « 0.425). Note finally that the substantial difference 
between the LIM and FCI excitation energies prevents 
the effective DD and shifted auxiliary excitation energy 
curves to be symmetric with respect to the weight axis, 
as it should be in the exact theory. 


£. Range-separated results for He 

As illustrated in the middle and bottom panels of 
Fig. 3, the auxiliary excitation energy, shown for ^ = 0.4 
and l.Oa]/^, becomes linear in ru as /i increases. This is 
in agreement with the first-order derivative expression 
in Eq. (B7). Indeed, when /i —> -l-oo, the auxiliary 
wavefunctions become the physical ones which are 
w-independent. Consequently, the third term on the 
right-hand side, that is responsible for the minimum 
at w = 0.01 observed when /i = 0, vanishes for larger 
pL values. Similarly, the auxiliary energies will become 
ru-independent and equal to the physical energies, thus 
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leading to a iw-independent first-order derivative. Inter¬ 
estingly, the (negative) second term on the right-hand 
side of Eq. (B7) is quadratic in the short-range kernel 
and is taken into account only when calculations are 
performed self-consistently. Since the short-range kernel 
becomes small as ^ increases, it is not large enough to 
compensate the positive contribution from the first term 
that is linear in the short-range kernel. As a result, the 
slope of the auxiliary excitation energy is positive for all 
w values. It also becomes clear that self-consistency will 
decrease the slope. 

Turning to the LIM excitation energies and the 
effective DDs, the former become closer to the FCI value 
as /i increases while the latter are reduced, as expected. 
The fact that the auxiliary excitation energy equals the 
LIM one for w = 0.25 confirms that the range-separated 
ensemble energy is essentially quadratic in w when 
^ Even though no accurate values for the 

short-range DD are available in the literature for any ru. 
Fig. 2 in Ref. [37] provides reference values for w = 0 
that are about 0.008 and 0.005 for /r = 0.4 and 
I.Ouq respectively. These values are simply obtained 
by subtracting the auxiliary excitation energies (denoted 
in Ref. [37]) from the standard FCI value (^ —)■ -l-oo 
limit). The effective DDs computed at the srLDA level 
for /i = 0.4 and I.Ooq ^ differ from these reference values 
by about a factor of ten. Note that srLDA and srPBE 
functionals give very similar results. 


3. Be and the stretched HeH'^ moleeule 

GOK-LDA and srLDA (/i = 0.4 and I.OUfC^) results 
are presented for Be and the stretched HeH+ molecule 
in Fig. 4. In both systems, the ensemble contains the 
ground state and a first singly-excited state, exactly like 
for He. Effective DD curves share similar patterns but 
their interpretations differ substantially. Let us first 
consider the Be atom. At the GOK-LDA level (top 
left-hand panel in Fig. 4), self-consistency effects are 
important. They are responsible for the negative slope of 
the auxiliary excitation energy at w = 0. Interestingly, 
the slope at w = 0 is larger in absolute value for He 
than for Be. This is clearly shown in the bottom panel 
of Fig. 5. As the auxiliary excitation energy decreases 
on a broader interval than for He, the second term on 
the right-hand side of Eq. (BIO) might become larger in 
absolute value as w increases. Its combination with the 
third term (linear in w) may explain why the minimum 
is reached at a larger ensemble weight value than for 
He {w « 0.045). One may also argue that this third 
term, that is only described at the self-consistent level, 
is smaller for Be than for He, thus leading to a less 
pronounced curvature in w, as shown in the top panel 
of Fig. 5. The auxiliary excitation energy becomes 
linear in w when /x = 0.4 and I.Ooq ^ (see middle and 


bottom left-hand panels in Fig. 4). Note finally that 
the effective DDs are about ten times smaller than in He. 

Let us now focus on the stretched HeH+ molecule. As 
shown in Fig. 5, patterns observed at the GOK-LDA level 
for He and Be are strongly enhanced due to the charge 
transfer. The interpretation is however quite different. 
Indeed, as shown in the top right-hand panel of Fig. 4, 
self-consistency is negligible for small w values and is 
therefore not responsible for the large negative slope of 
the auxiliary excitation energy at w = 0. This was ex¬ 
pected since the self-consistent contribution to the slope 
(second term on the right-hand side of Eq. (BIO)) in¬ 
volves the overlap between the HOMO (localized on He) 
and the LUMO which is, in this particular case, strictly 
zero. Consequently, as readily seen in Eq. (BI2), the 
(negative) LDA exchange and correlation kernels [3] are 
responsible for the negative slope at w = 0. The lat¬ 
ter is actually smaller in absolute value when the LDA 
correlation density functional is set to zero in the cal¬ 
culation (not shown), thus confirming the importance of 
both exchange and correlation contributions to the ker¬ 
nel. Note that, as w increases, self-consistency effects are 
growing. This can be related with the third term on the 
right-hand side of Eq. (BIO) where the response of the ex¬ 
cited state to changes in w contributes. Interestingly, for 
/X = 0.4a(C^, the contribution to the slope, at ic = 0, from 
the short-range exchange-correlation kernel is significant 
enough [3] so that the pattern observed at the GOK- 
LDA level does not completely disappear (see the middle 
right-hand panel in Fig. 4). On the other hand, for the 
larger /i = I.Ouq ^ value, the auxiliary excitation energy 
becomes essentially linear in w with a positive slope (see 
the bottom right-hand panel in Fig. 4). Note finally that 
the stretched HeH+ molecule exhibits the largest effec¬ 
tive DDs. 


4- H2 

Results obtained for H 2 are shown in Figs. 5 and 6. 
At equilibrium, they are quite similar to those obtained 
for He. Still, at the GOK-LDA level, the negative 
slope of the auxiliary excitation energy at rc = 0 is not 
related with self-consistency (see the top left-hand panel 
in Fig. 6), in contrast to He. Self-consistency effects 
become significant as w increases. Effective DDs at 
w = 0 are equal to 40.9, 36.2 and 8.6 mEh for /x = 0, 0.4 
and l.Oa)]'^, respectively. They are significantly larger 
than the accurate values deduced from Fig. 6 in Ref. [37] 
(7.1, 5.7 and about zero vtiEh). 

In the stretched geometry (right-hand panels in Fig. 6), 
the nature of the first excited state completely changes. 
It corresponds to the double excitation —)■ Icr^. At 

the GOK-LDA level, self-consistency effects are negligi¬ 
ble. This was expected since, according to Eq. (B7), the 
latter effects involve couplings between ground and ex- 
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cited states through the density operator. Consequently, 
a doubly-excited state will not contribute. Moreover, 
the difference in densities between the ground-state and 
first doubly-excited GOK determinants reduces along the 
bond-breaking coordinate, simply because the overlap 
between the Is orbitals reduces. As a result, the first- 
order derivative of the auxiliary excitation energy is very 
small, as confirmed by Fig. 5. This analysis holds also 
for larger fj, values. The only difference is that, when 
/i > 0, both ground- and excited-state wavefunctions are 
multiconfigurational [54, 55]. In a minimal basis, they 
are simply written as 



In this case, both ground and excited states have the 
same density, 

(r) = (r) = ^ (n ^2 (r) -F n ^2 (r)) , (108) 

and their coupling through the density operator equals 

(^ol^(i’)l^'i) = ^(n„2(r) -n<,2(r)), (109) 

which is zero as the overlap between the Is orbitals is 
neglected. 

Since the ensemble energy is, for any fj, value, almost 
linear in w, the LIM and auxiliary excitation energies 
are very close for any weight. Consequently, the effective 
DD is very small (4.5 mEh for fj, = Ooq^ and w = 0). 
Since the deviation of the LIM excitation energy from 
the FCI one is relatively large (about —0.12Efi for fj, = 
Ooq ^), symmetry of the plotted curves with respect to 
the weight axis is completely broken, in contrast to the 
other systems. In this particular situation, LIM brings 
no improvement and the effective DD is expected to be 
far from the true DD. For comparison, the latter equals 
about 200 mEh for a slightly larger bond distance (4.2ao) 
and ^ = Ooq according to Fig. 7 in Ref. [37]. For the 
same distance, the KS-LDA auxiliary excitation energy 
(not shown) deviates by ISOmAft, in absolute value from 
the FCI value, which is in the same order of magnitude 
as the true DD. Therefore, for R = 3.7ao, the true DD is 
expected to be much larger than the effective one. 

B. Excitation energies 

1. Single excitations 

LIM excitation energies have been computed when 
varying /x for the various systems studied previously. Sin¬ 
gle excitations are discussed in this section. Results are 
shown in Fig. 7. It is quite remarkable that, already 
for fi = 0, LIM performs better than standard TD-DFT 


with the same functional (LDA or PBE). This is also true 
for the 2S+ charge transfer state in the stretched HeH+ 
molecule. We even obtain slightly better results than 
the popular TD-CAM-B3LYP method. As expected, the 
error with respect to FCI reduces as /i increases. Note 
that, for He, it becomes zero and then changes sign in 
the vicinity of ^ = LOug The latter value gives also 
accurate results for the other systems, which is in agree¬ 
ment with Ref. [13]. Note also that, for the typical value 
^ = 0.4 — 0.5ag ^ [25, 26], the slope in /x for the LIM 
excitation energy is quite significant. It would there¬ 
fore be relevant to adapt the extrapolation scheme of 
Savin [36, 38] to range-separated ensemble DFT. This 
is left for future work. Note that srLDA and srPBE 
functionals give rather similar results. For comparison, 
auxiliary excitation energies obtained from the ground- 
state density {w = 0) are also shown. The former re¬ 
duce to KS orbital energy differences for /x = 0. In this 
case, TD-DFT gives slightly better results, except for the 
charge transfer excitation in HeH+ where the difference 
is very small, as expected [1]. Both srLDA and srPBE 
auxiliary excitation energies reach a minimum at rela¬ 
tively small /X values (0.125ag ^ for He). This is due to 
the approximate short-range (semi-) local potentials that 
we used. Indeed, as shown in Ref. [37], variations in /x 
are expected to be monotonic for He and H 2 at equi¬ 
librium if an accurate short-range potential were used. 
Since the range-separated ensemble energy can be ex¬ 
pressed in terms of the auxiliary energies (see Eq. (47)), 
it is not surprising to recover such minima for some LIM 
excitation energies. Let us finally note that the auxil¬ 
iary excitation energy converges more rapidly than the 
LIM one to the FCI value when /x increases from l.Oug 
For Be, convergences are very similar. As already men¬ 
tioned, the convergence can actually be further improved 
by means of extrapolation techniques [36, 38]. In conclu¬ 
sion, the LIM approach is promising at both GOK-DFT 
and range-separated ensemble DFT levels. In the latter 
case, /X should not be too large otherwise the use of an 
ensemble is less relevant. Indeed, auxiliary excitation en¬ 
ergies obtained from the ground-state density are in fact 
better approximations to the FCI excitation energies, at 
least for the systems and approximate short-range func¬ 
tionals considered in this work. This should be tested on 
more systems in the future. 


2. Double excitations 

One important feature of both GOK and range- 
separated ensemble DFT is the possibility of modeling 
multiple excitations, in contrast to standard TD-DFT. 
Results obtained for the and l^D states in 

the stretched H 2 molecule and Be, respectively, are 
shown in Fig. 8. We focus on H 2 first. As discussed 
previously, LIM and auxiliary excitation energies are 
almost identical in this case. For = Oog they differ 
by about -0.12 Eh from the FCI value. There are 
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no significant differences between srLDA and srPBE 
results. The error monotonically reduces with increasing 
fi. Interestingly, for ^ = O.doQ^, the LIM excitation 
energy equals 0.2i7Eh, that is very similar to the 
multi-configuration range-separated TD-DFT result 
obtained with the same functionals (0.238Eh). [4] This 
confirms that the short-range kernel does not contribute 
significantly to the excitation energy, since the ground 
and doubly-excited states are not coupled by the density 
operator (see Eq. (109)). Note that, for R = 4.2ao 
and fj, = 0.4a)C^, the srLDA auxiliary excitation energy 
(not shown) equals 0.194E?i, that is rather close to 
the accurate value (O.lSlE/j) deduced from Fig. 7 in 
Ref. [37]. As a result, the approximate (semi-)local 
density-functional potentials are not responsible for the 
large error on the excitation energy. One would blame 
the adiabatic approximation if TD linear response theory 
were used. In our case, it is related to the WIDFA 
approach. In this respect, it seems essential to develop 
weight-dependent exchange-correlation functionals for 
ensembles. Applying the GAGE formalism to model 
systems would be instructive in that respect. Work is 
currently in progress in this direction. 

Turning to the doubly-excited l^D state in Be, LIM 
is quite accurate already at the GOK-DFT level. In¬ 
terestingly, the largest and relatively small errors in ab¬ 
solute value (about 4.0 and 7.0 uiEh for the srLDA 
and srPBE functionals, respectively) are obtained around 
/i = l.Oa])'^. In this case, the ensemble contains four 
states (1^5', 2^S and two degenerate 1^1? states) whereas 
in all previous cases first excitation energies were com¬ 
puted with only two states. This indicates that /i values 
that are optimal in terms of accuracy may depend on 
the choice of the ensemble. This should be investigated 
further in the future. 


V. CONCLUSIONS 

A rigorous combination of wavefunction theory with 
ensemble DFT for excited states has been investigated 
by means of range separation. As illustrated for simple 
two- and four-electron systems, using local or semi-local 
ground-state density-functional approximations for 
modeling the short-range exchange-correlation energy 
of a bi-ensemble with weight w usually leads to range- 
separated ensemble energies that are not strictly linear 
in w. Consequently, the approximate excitation energy, 
that is defined as the derivative of the ensemble energy 
with respect to w, becomes w-dependent, unlike the 
exact derivative. Moreover, the variation in w can be 
very sensitive to the self-consistency effects that are 
induced by the short-range density-functional potential. 

In order to define unambiguously approximate ex¬ 
citation energies in this context, we proposed a linear 
interpolation method (LIM) that simply interpolates 


the ensemble energy between w = 0 (ground state) and 
w = 1/2 (equiensemble consisting of the non-degenerate 
ground and first excited states). A generalization to 
higher excitations with degenerate ground and excited 
states has been formulated and tested. It simply consists 
in interpolating the ensemble energy linearly between 
equiensembles. LIM is applicable to GOK-DFT that is 
recovered when the range-separation parameter /r equals 
zero. In the latter case, LIM performs systematically 
better than standard TD-DFT with the same functional, 
even for the 2S+ charge-transfer state in the stretched 
HeH+ molecule. For typical values ^ = 0.4 — 0.5aQ 
LIM gives a better approximation to the excitation 
energy than the auxiliary long-range-interacting one 
obtained from the ground-state density. However, for 
larger fj, values, the latter excitation energy usually 
converges faster than the LIM one to the physical result. 

One of the motivation for using ensembles is the 
possibility, in contrast to standard TD-DFT, to model 
double excitations. Results obtained with LIM for the 
1^1? state in Be are relatively accurate, especially at the 
GOK-DFT level. In the particular case of the stretched 
H 2 molecule, the range-separated ensemble energy is 
almost linear in w, thus making the approximate 
excitation energy almost weight-independent. LIM 
brings no improvement in that case and the error on 
the excitation energy is quite significant. This example 
illustrates the need for weight-dependent exchange- 
correlation functionals. Combining adiabatic connection 
formalisms [14] with accurate reference data [15] will 
hopefully enable the development of density-functional 
approximations for ensembles in the near future. 

Finally, in order to turn LIM into a useful mod¬ 
elling tool, a state-averaged complete active space self- 
consistent field (SA-CASSCF) should be used rather than 
Cl for the computation of long-range correlation effects. 
Since the long-range interaction has no singularity at 
ri 2 = 0, we expect a limited number of configurations 
to be sufficient for recovering most of the long-range cor¬ 
relation. This observation has already been made for the 
ground state [33, 56]. Obviously, the active space should 
be chosen carefully in order to preserve size consistency. 
The implementation and calibration of such a method is 
left for future work. 
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Appendix A: SELF-CONSISTENT 
RANGE-SEPARATED ENSEMBLE 
DENSITY-FUNCTIONAL PERTURBATION 
THEORY 

The self-consistent Eq. (45) can be solved for small 
w values within perturbation theory. For that purpose 
we partition the long-range interacting density-functional 
Hamiltonian as follows, 

(Al) 

where, according to Eq. (9), the perturbation equals 


Consequently, 

dn^ = (1 - - n°) 

+ 00 

= ^.F'=.F(n>,.-nO) 

= - n°) + ... (A9) 


Appendix B: DERIVATIVE OF THE AUXILIARY 
EXCITATION ENERGY 


= / dr| ]^(r),(A2) 


(5n(r) 

and, according to Eq. (44), 


i5n(r) 


h^.“(r) =n°(r)-b w 




dw 


+ 0{w^) 


\ w—0 

n° (r) -I- w (r) — nP (r 

an^M.»(r) 


-b W- 


dw 


+ 0{w^). (A3) 




Combining Eq. (A2) with Eq. (A3) leads to 

-b o(w) 


drdr'-^— 

5n{v')5n{r 


y(^n,j,M(r') -n°(r') 


dnq,i.^{v') 


dw 


h(r) + 0{w). (A4) 


I LP = 0 ' 

From the usual first-order wavefunction correction ex¬ 
pression 




dw 


u=0 j>l 


si;-sit ■ 


(A5) 


and the expression for the derivative of the ground-state 
density, that we simply denote dn^, 


dn^{v^) = 


9n^f;.™(ri) 


dw 


w—Q 


= 2(1'; 




dw 


10=0 


we obtain the self-consistent equation 

dn^ = Pdn^ -b P (n,i,M — n°), 


(A6) 


(A7) 


where is a linear operator that acts on any function 
/(r) as follows. 


ff(w = 2Y, 


i>i 


6n{r')6n{r) S^ — St 


According to Eq. (48), the first-order derivative of the 
individual auxiliary energies can be expressed as 


dSt 


dw 


dr'dr i 

,5n(rO<5n(r) 

ah^’“(r') 


dw 


n^^.™(r), (Bl) 


where 


dw 


= 5h^’’"(r')-b 


an^M,™(r') 


dw 


-bW 


dSn^’^{r') 
dw ' 


and 


(5h^’'"(r') = n^M.»(r') - n^^.™(r'), 


(B2) 


(B3) 


so that the derivative of the auxiliary excitation energy 
in Eq. (49) can be written as 


dAf'^’' 

dw 


dn{r')6n{r) 


X j^(5h''’“(r')(5h'^’“(r) + 

9(5h^’“(r') 


an^M,™(r') 


dw 




+w 


dw 


-5h^’“'(r) . 


(B4) 


According to perturbation theory through first order (see 
Appendix A), the response of the ground-state density to 
variations in the ensemble weight equals 


(9n^^,™(r') 

dw 


= 2( 4-^’' 


h(r') 




= 2 E 


i>i 


dridr2 


dw 


i5n(r2)i5n(ri) 


UK_ 


r(r')<r(ri)5hyyy(r2) 

dw 


, (B5) 


<2(1’) = (^ol^WI^D- 


(AS) 


where ng-“’(r') = (^g’“'|h(r')|'Ff’“). Note that, as al¬ 
ready pointed out for ic = 0 (see Eq. (A7)), Eq. (B5) 
should be solved self-consistently. By considering the 
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first contribution to the response of the ensemble den¬ 
sity in Eq. (B2) we obtain 


dw 


= 2 E 


i>l 


dridr2 


5n{r-2)5n{Yx) 




Sh^’'"ir2) 


(B6) 


thus leading to the following expansion 


dw 

+ 2 E 


r2 

dr'dr ^ , E fi^’’^(r')(5h^’^(r) 

(3n(r')dn(r) 


1 


nl-L,w cMi'^ 


dr'dr , (rXr (r') 


-l-w 


5n{r')5n{v) 


dr dr 


dn{r')6n{r) dw 


,5h^’“(r) 


+ ... (B7) 

Note that, at the srLDA level of approximation, the 
exchange-correlation contribution to the short-range ker¬ 
nel is strictly local [3]. By using the decomposition 

X2 psrLDA.^r 1 


6n{r')Sn{r) 


dri^ 


Sir-r'), (B8) 


the first term on the right-hand side of Eq. (B7) can be 
simplihed as follows, 

c2 ^srLDA,Air~/i,i«i 

dr'dr- XXEa— 

dn{r)dn[r) 

J J dr'drwX(|i--r'|)(5h'^’“(r')(5h^’’"(r) 


/ an2 


((5h''’“(r))' 


(B9) 


In the GOK-DFT limit (/i = 0), if the first excitation 
is a single excitation from the HOMO to the LUMO, 
the auxiliary excitation energy reduces to an orbital en¬ 
ergy difference Ae™ whose derivative can formally be ex¬ 
pressed as follows, according to Eq. (B7), 


dAe*^ 

dw 

-k4 


1 / I ^ ^Hxc[fi ]c~w/ l\X~w/ \ 

dr dr——————dn (r )6n (r) 

on(r )dn(r) 


E 


i<N/2,a>N/2 


e- — £„ 


+W 


dn{r)dn[r) 

drdr—————— - Sn r) 

dn{r)dn{r} ow j 


(BIO) 


where 


A/2-1 

h“(r) = 2 ^ C(r)^ 

k^l 

+ (2 - w)X/2W^ + W^N/2+li^f^ 

6n'"{r) = X/ 2 +iW^ - 


and {0“(i’)}fe are the GOK-DFT orbitals with the associ¬ 
ated energies {s^}k that are obtained within the WIDFA 
approximation. Note that, in practical calculations, par¬ 
tially occupied GOK-DFT orbitals have not been com¬ 
puted explicitly. Instead, we performed FGI calculations 
in the basis of determinants constructed from the KS or¬ 
bitals. 

Let us hnally note that if the HOMO and LUMO do not 
overlap, the first term on the right-hand side of Eq. (BIO) 
can be further simplihed at the LDA level, according to 
Eq. (B9), thus leading to 


drdr , , P. , . Sn (r)6n (r) 
dn{r)dn(r) 

, , , X/2(r)"X/2(l’')" 

dr dr —^/- 


r — r| 

2lw t^r\2 


dr 


, , , X/2+l(l’)"</>A/2+l(r') 

dr dr — - -;- y. - 

|r - r'l 

^^exc(fd"(r)) 

d'n? 


x(</'a/2W +^N/2+i{^Y 


(B12) 


Appendix C: SELF-CONSISTENCY EFFECTS ON 
THE ENSEMBLE AND AUXILIARY ENERGIES 

Let n denote a trial ensemble density for which the 
auxiliary wavefunctions can be determined: 

XhixN) = CWIXW), * = o,i. (Cl) 

The resulting auxiliary ensemble density, 

n“[n](r) = (1 - w)nv[,M[„](r)-h wnvi/ 5 ‘[„](r), (^2) 

is then a functional of n, like the ensemble energy that 
can be expressed as 

E^’'“[n] = (1 — w)£1q [n] -I- w5f [n] 

/ r r 1 

^ W]- (C3) 

The converged ensemble density fulhls the following 
condition: 


(C4) 

If we now consider variations around the trial density, 
n ^ n + 6n, the ensemble energy will vary through hrst 


+ ... 
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order in 6n as follows, 


SE^''^[n] = (1 — w)6EQ[n] + w6Ei[n] 


./dr^( 


SE\ 


isr,/j. r 
Hxc i 


dr 


y 5n(r) 

6E^^^[n]] 

Sn{r) 



(C5) 


where, according to the Hellmann-Feynman theorem. 



Combining Eqs. (Cl) and (C5) with Eq. (C6) leads to 





Sn{r) 


SE^\ 

Sn{r) j 


X Jn™ [n] (r). (C7) 


According to Eq. (C6), the auxiliary excitation energy 
AS^ [n] = Si [n] — Sq [n] will vary through first order as 


6AS^^[n] = 


drdr' Sniv') 

5n{v')5n{vy ^ ^ 


(C8) 


We conclude from Eqs. (C4), (C7) and (C8) that vari¬ 
ations 6n around the converged ensemble density 
will induce at least first and second order deviations in 
Sn for the auxiliary excitation and ensemble energies, re¬ 
spectively. 
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FIGURE CAPTIONS 

Figure 1: (Color online) Range-separated ensemble en¬ 
ergy obtained for He at the WIDFA level when 
varying the ensemble weight w for fj, = 0 and 
l.OaQ^. Comparison is made with the linear in¬ 
terpolation method (LIM) for /x = and FCI. 

The ensemble contains both and 2^S states. 
The srLDA functional has been used. 

Figure 2: (Color online) Schematic representation of 
the linear interpolation method. Ensemble ener¬ 
gies and their first-order derivatives are shown in 
the top and bottom panels, respectively. See text 
for further details. 

Figure 3: (Color online) Effective DD (A^^"), auxiliary 
(A£i^’“) and LIM (w^j,^) excitation energies asso¬ 
ciated with the excitation 1^5 —)■ 2^5” in He. Re¬ 
sults are shown for fj, = 0, 0.4 and 1.0 with 
the srLDA (left-hand panels) and srPBE (right- 
hand panels) functionals when varying the ensem¬ 
ble weight w. Comparison is made with the FCI ex¬ 
citation energy wpci = 0.7668 Eh- Empty squares 
are used for showing non-self-consistent results. 

Figure 4: (Color online) Effective DD (A^^"), auxiliary 
(A£i^’“) and LIM (w^j,^) excitation energies asso¬ 
ciated with the excitations 1^5 —)■ 2^5' in Be (left- 
hand panels) and 1^S+ —)■ 2^E+ in the stretched 
HeH+ molecule (right-hand panels). Results are 
shown for = 0, 0.4 and I.Ooq ^ with the srLDA 
functional when varying the ensemble weight w. 
Comparison is made with the FCI excitation ener¬ 
gies (wpci = 0.2487i?;i for Be and wpci = 0.4024if^ 
for HeH+). Empty squares are used for showing 
non-self-consistent results. 

Figure 5: (Color online) Auxiliary excitation energies 
obtained with /i = Ouq ^ and the srLDA functional 
(that is equivalent to GOK-LDA) when varying the 
ensemble weight w in the various systems consid¬ 
ered in this work. See text for further details. Exci¬ 
tation energies are shifted by their values at w = 0 
for ease of comparison. A zoom is made on the 
0 < ru < O.I region in the bottom panel. 

Figure 6: (Color online) Effective DD (A((g"), auxiliary 
(A5^’“) and LIM (w^j,^) excitation energies asso¬ 
ciated with the excitation —)■ 2^E+ in H 2 at 

equilibrium (left-hand panels) and in the stretched 
geometry (right-hand panels). Results are shown 
for /i = 0, 0.4 and I.0a(C^ with the srLDA functional 
when varying the ensemble weight w. Comparison 
is made with the FCI excitation energies (wpci = 
0.4828if;i at equilibrium and wpci = 0.3198if?i in 
the stretched geometry). Empty squares are used 
for showing non-self-consistent results. 


Figure 7: (Color online) LIM excitation energies ob¬ 
tained for the single excitations discussed in this 
work with srLDA and srPBE functionals when 
varying the range-separation parameter /x. Com¬ 
parison is made with standard TD-DFT and FCI. 
For analysis purposes, auxiliary excitation energies 
obtained from the ground-state density {w = 0) are 
shown (curves with empty circles). 

Figure 8 : (Color online) LIM excitation energies cal¬ 
culated for the doubly-excited 2^E+ state in the 
stretched H 2 molecule (top panel) and 1^1? state 
in Be (bottom panel) when varying the range- 
separation parameter /x with srLDA and srPBE 
functionals. Comparison is made with FCI. For 
H 2 , auxiliary excitation energies obtained from the 
ground-state density {w = 0) are shown (curves 
with empty circles) for comparison. 
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FIG. 2: Senjean et al, Phys. Rev. A 
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FIG. 3: Senjean et al, Phys. Rev. A 
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FIG. 4: Senjean et al, Phys. Rev. A 
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FIG. 5: Senjean et al, Phys. Rev. A 
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FIG. 6: Senjean et al, Phys. Rev. A 
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FIG. 7: Senjean et al, Phys. Rev. A 
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FIG. 8: Senjean et al, Phys. Rev. A 
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